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ABSTRACT 


This  paper  reports  on  an  experiment  In  trying  to  understand  an 
unfamiliar  program  of  some  complexity  and  to  record  the  authors'  under¬ 
standing  of  It.  The  goal  was  to  simulate  a  practicing  programmer  in  a 
program  maintenance  environment  using  the  techniques  of  program  design 
adapted  to  program  understanding  and  documentation;  that  is,  given  a 
program,  a  specification  and  correctness  proof  were  developed  for  the 
program.  The  approach  points  out  the  value  of  correctness  proof  ideas  in 
guiding  the  discovery  process.  Toward  this  end,  a  variety  of  techniques 
were  used:  direct  cognition  for  smaller  parts,  discovering  and  verifying 
loop  Invariants  for  larger  program  parts,  and  functions  determined  by 
additional  analysis  for  larger  program  parts.  An  indeterminate  bounded  variable 
was  introduced  into  the  program  documentation  to  summarize  the  effect  of 
several  program  variables  and  simplify  the  proof  of  correctness. 
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UNDERSTANDING  AND  DOCUMENTING  PROGRAMS 


I.  INTRODUCTION 

Understanding  programs  -  We  report  here  on  an  experiment  in  trying  to 
understand  an  unfamiliar  program  of  some  complexity  and  to  record  our  under¬ 
standing  of  It.  We  are  as  much  concerned  with  recording  our  understanding  as 
with  understanding.  Every  day  programmers  are  figuring  out  what  existing 
programs  do  more  or  less  accurately.  But  most  of  this  effort  is  lost,  and 
repeated  over  and  over,  because  of  the  difficulty  of  capturing  this  under¬ 
standing  on  paper.  We  want  to  demonstrate  that  the  very  techniques  of  good 
program  design  can  be  adapted  to  problems  of  recording  hard  won  understandings 
about  existing  programs. 

In  program  design,  we  advocate  the  joint  development  of  design  and  correct¬ 
ness  proof,  as  shown  by  Dljkstra  in  (Dahl,  Dljkstra,  and  Hoare)  and  (Dijkstra) 
and  by  (Linger,  Mills,  and  Witt),  rather  than  a  posteriori  proof  development. 
Nevertheless,  we  believe  that  the  idea  of  program  correctness  provides  a  com¬ 
prehensive  a  posteriori  strategy  for  developing  and  recording  an  understanding 
of  an  existing  program.  In  fact,  we  advocate  another  kind  of  joint  develop¬ 
ment,  this  time,  of  specification  and  correctness  proof.  In  this  way,  we  have 
a  consistent  approach  dealing  always  with  three  objects;  namely,  (1)  a  specifi¬ 
cation,  x2)  a  program,  and  (3)  a  correctness  proof.  In  writing  a  program,  we 
are  given  (1)  and  develop  (2)  and  (3)  jointly;  in  reading  a  program,  we  are 
given  (2)  and  develop  (1)  and  (3)  jointly.  In  either  case,  we  end  up  with  the 
same  harmonious  arrangement  of  (1)  and  (2)  connected  by  (3)  which  contains  our 
understanding  of  the  program. 

In  the  experiment  at  hand,  our  final  understanding  exceeded  our  most 
optimistic  initial  expectations,  even  though  we  have  seen  these  ideas  succeed 
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before.  One  new  Insight  from  this  experiment  was  how  little  we  really  had  to 
know  about  the  program  to  develop  a  complete  understanding  and  proof  of  what 
it  does  (in  contrast  to  how  it  does  it).  Without  the  correctness  proof  ideas 
to  guide  us,  we  simply  would  not  have  discovered  how  little  we  had  to  know. 

In  fact,  we  know  a  great  deal  acre  than  we  have  recorded  here  about  how  the 
program  works,  which  we  chalk  up  to  the  usual  dead  ends  of  a  difficult  discovery 
process.  But  the  point  is,  without  the  focus  of  a  correctness  proof,  we  would 
still  be  trying  to  understand  and  record  a  much  larger  set  of  logical  facts 
about  the  program  than  is  necessary  to  understand  precisely  what  it  does. 

In  retrospect,  we  used  a  variety  of  discovery  techniques.  For  simpler 
parts  of  the  program,  we  used  direct  cognition.  In  small  complex  looping  parts, 
we  discovered  and  verified  loop  Invariants.  In  the  large,  we  organized  the 
effect  of  major  program  parts  as  functions  to  be  determined  by  additional 
analysis.  We  also  discovered  a  new  way  to  express  the  effect  of  a  complex 

program  part  by  introducing  a  bounded  Indeterminate  variable  which  radically. _ 

simplified  the  proof  of  correctness  of  the  program  part. 

The  experiment  -  We  were  interested  in  a  short  but  complex  program  using 
real  arithmetic,  and  felt  that  more  attention  might  be  paid  to  the  structure 
and  correctness  of  programs  that  deal  with  real  arithmetic.  The  program  was 
chosen  by  Professor  James  Vandergraft  of  the  University  of  Maryland  as  a  diffi¬ 
cult  program  to  understand.  It  was  a  FORTRAN  program  called  ZEROIN  which 
claimed  to  find  a  zero  of  a  function  given  by  a  FORTRAN  subroutine. 

Our  goal  was  to  simulate  a  practicing  programmer  in  a  program  maintenance 
environment.  We  were  given  the  program  and  told  its  general  function.  The 
problem  then  was  to  understand  it,  verify  its  correctness,  and  possibly 
modify  it,  to  make  it  more  efficient  or  extend  its  applicability.  We  were  not 
given  any  more  about  the  program  than  the  program  itself.  The  program  given 


Co  us  la  shown  in  Figure  1.  Professor  Vander graft  played  Che  role  of  a  user 
of  Che  program  and  posed  four  questions  regarding  the  program: 

1.  I  have  a  lot  of  equations,  some  of  which  might  be  linear.  Should 
I  test  for  linearity  and  then  solve  the  equation  directly,  or  just 
call  ZEROIN?  That  is,  how  much  work  does  ZERO IN  do  to  find  a  root 
of  a  linear  function? 

2.  What  will  happen  if  I  call  ZEROIN  with  FA  and  FB  both  positive? 

How  should  the  code  be  changed  to  test  for  this  condition? 

3.  It  Is  claimed  that  the  inverse  quadratic  Interpolation  saves  only 

.5  function  evaluations  on  the  average.  To  get  a  shorter  program,  I 
would  like  to  remove  the  inverse  quadratic  interpolation  part  of  the 
code.  Can  this  be  done  easily?  How? 

4.  Will  ZEROIN  find  a  triple  root? 

It  should  be  noted  that  the  authors  are  not  currently  working  in  the  area  of 
numerical  analysis,  though  it  is  not  an  unknown  area  to  them. 


ZEROIN. PROGRAM 
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REAL  FUNCTION  ZEROIN  (AX,  3X,  F,  TCL,IP) 

RcAL  AX,  3 X ,  F,  TOL 

REAL  A,  3,  C,  0,E,  EPS,  FA,  F3,  FC ,  TOLl, 

COM’JTS  EPS,  THE  RELATIVE  MACHINE  »RECISION 

EPS  =  1.3 
13  EPS  =  EPS/2.2 
mi  *  1.3  ♦  EPS 
IF  (TOLl  .GT.  1.3  )  GO  TO  ID 

I M I T I ALIZATION 

IF  (IP  .EG.  1)  waiTE(6,11> 

FORMAT('  THE  INTERVALS  DETERMINED  3Y  ZEROIN  ARE") 
A  *  AX 
3  =  3X 
FA  a  F ( A ) 

F 3  =  F  (3) 

BEGIN  STEP 


23  C  =  A 
FC  =  FA 
3  =  a  -  a 
E  *  D 

33  IF  (IP  .EG.  1)  WRITE(5,21>  3,C 

31  FORMAT  (2E15.3) 

IF  (  A53(FC)  . 3  E .  A3  S ( F3 )  )  GO  TO  43 
A  a  3 
3  a  c 

C  a  A 
FA  a  pa 
F  3  *  fc 
FC  *  FA 


C0N/ER3ENCE  TEST 

43  TOLl  *  2 .3*EPS*A3S (3)  ♦  2.5*T0L 
<M  a  . 5*  ( c  -  3) 

IF  (A3S(XM)  .LS.  TOLl  )  30  TO  / D 


IF  (  F3 


S.3  )  30  TO  9 C 


IS  BISECTION  NECESSARY 

IF  (  A  3  S  (  E  )  .LT.  TOLD  GO  TO  73 
IF(A3S(FA)  .LE.  A3S(F3)  )  30  TO  7G 

IS  1JADRATIC  INTERPOLATION  POSSIBLE 

IF  (A  .NE.  C)  30  TO  5  C 

LINEAR  INTERPOLATION 

3  =  F3/FA 
»  *  2.3*XM*S 
1  »  1.3  -  S 
30  TO  53 

INVERSE  JUADRATi:  INTERPOLATION 

53  J  *  FA/FC 
R  *  F3/FC 
3  *  r  3  /F  A 

»  a  S*(2.D*X.M*3*(3  -  R)  -  (3  -  A )  *  ( R  -  1.2)) 
3  *  (3  -  1  .D)*(R  -  «  ,3)*(S  -  1 ,D) 

A3JJST  SIGNS 


53  IF  (P  .GT.  3.3  ) 
»  *  A?  S  (  P  ) 
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9  7  . 
93  . 
9?  . 
133  . 
101  . 
132  . 


C 

C  IS  INTERPOLATION  ACCS  PT  A3  LE 

C  If  <(2.3*P>  .3=.  <3.C*XM*S  -  A3S(T0L1*3>>>  SO  TO  ?3 

IF  (P  «3£  .  A3S(3.5*E*Q)>  GO  TO  73 
S  *  0 
o  =  p/a 
30  TO  33 
C 

C  BISECTION 
C 

73  0  «  X* 

3  =  D 
C 

C  COM»LETE  STEP 
C 

S3  A  =  3 
FA  =  F  3 

IF  ( A  3  S ( 0 )  .GT.  T0L1  )  £  =  3  ♦  9 

IF  (  A  3  S  (  D  )  .LS.  TOLD  £  =  3  ♦  S I S N ( T OL 1 , XM ) 

FB  =  FC3) 

IF  <  (F3* (FC/A3S(FC ) )  3  .3T.  3.3  >  GO  TO  20 
30  TO  33 

V* 

C  DONS 
C 

e3  ZEROIN  *  3 
RETURN 
END 


*****  ZEROIN. INFO  ***** 


ZEROIN  IS  A  FUNCTION  SUBPROGRAM  WHICH  FINDS 
A  ZERO  OF  THE  FUNCTION  F(X)  IN  THE  INTERVAL  AX,?X  . 

THE  CALLING  STATEMENT  SHOULD  HAVE  HE  FOR* 

X*  =  ZER0IN(AX,3X,F,TOL,IP> 

•HERE  THE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS. 

INPJT 

AX  LEFT  ENDPOINT  OF  INITIAL  INTERVAL 
3  <  RI3HT  ENDPOINT  OF  INITIAL  INTERVAL 

F  FUNCTION  SUBPROGRAM  WHICH  EVALUATES  F(X)  FOR  ANT  X  IN 

THE  INTERVAL  A  X  »  3  X 

TOL  DESIRED  LENGTH  OF  THE  INTERVAL  OF  UNCERTAINTY  OF  THfc 
FINAL  RESULT  (  ,3E.  2.3) 

1 9  AN  INTEGER  PRINT  FLAG.  WHEN  SET  TO  0,  NO  “RINTIVS 

V  ILL  3c  DONE  3  Y  ZEROIN.  IF  SET  TO  1,  THEN 
ALL  OF  THE  INTERVALS  COMPUTED  3T  ZEROIN  WILL 
3E  PRINTED  OJT  . 


OUTPUT 

ZEROIN  A3CISS*  APPROXIMATING  A  ZERD  OF  F  IN  THE  INTERVAL  AX.BX 


IT  IS  ASSUMED  THAT  F(AX)  AND  F(5X)  HAVE  OPPOSITE  S 1 3 N S 
WITHOJT  A  CHECK.  ZEROIN  RETURNS  A  ZERO  X  IN  THE  GIVEN  INTERVAL 
A  X  ,  3  <  TO  WITHIN  A  TOLERANCE  4 *M A C H E » S * A  3 S  (  X )  *  TOL,  WHERE  MACHEPS 
IS  THE  RELATIVE  MACHINE  PRECISION.  .  .. 

THIS  FUNCTION  SJ3“R0GRA,<!  IS  *  S^T  3HTLY  MODIFIED  TRANSLATION  OF 
THE  ALGOL  t3  PROCEDURE  ZERO  GIVEN  IN  RICHARD  SR  eNT  ,  AL3  OR  I THMS  FO1 
MINIMIZATION  WITHOUT  DERIVATIVES  ,  PRENTICE  -  HALL, INC.  (1C7:D. 
THIS  VERSION  IS  COPISO  F  R  ON  "COMPUTER  METHODS  FOR  MATHEMATICAL 
COMPUTATIONS'*  3Y  FORSYTHE,  NALCOLM,  AND  MOLEO.  THE  ONLY  CHANGE 
IS  THE  INCLUSION  OF  THE  PRINT  F L A 3  IP. 
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4 


II.  TECHNIQUES  FOR  UNDERSTANDING  PROGRAMS 

Flowcharts  -  Any  flowchartable  program  can  be  analyzed  in  a  way  we 
describe  next  £or  better  understandablllty  and  documentation.  For  a  fuller 
discussion,  see  (Linger,  Mills  and  Witt).  We  consider  flowcharts  as  directed 
graphs  with  nodes  and  lines.  The  lines  denote  flow  of  control  and  the  nodes 
denote  tests  and  operations  on  data.  Without  loss  of  generality,  we  consider 
flowcharts  with  just  three  types  of  nodes,  namely: 

function  node:  -4Th 

predicate  node: 

collecting  nodes: 

where  f  is  any  function  mapping  the  data  known  to  the  program  to  new  data, 
e.g.,  a  simple  FORTRAN  assignment  statement,  and  p  is  any  predicate  on  the  data 
known  to  the  program,  e.g.,  a  simple  FORTRAN  test.  An  entry  line  of  a  flowchart 
program  is  a  line  adjacent  to  only  one  node,  at  its  head;  an  exit  line  is 
adjacent  to  only  one  node,  its  tail. 

Functions  and  data  assignments  -  Any  function  mapping  the  data  known  to 
a  program  to  new  data  can  be  defined  in  a  convenient  way  by  generalized  forms 
of  data  assignment  statements.  For  example,  an  assignment,  denoted 
x  :-  e,  (e.g. ,  x  x  +  y) 

where  x  is  a  variable  known  to  the  program  and  e  is  an  expression  in  variables 
known  to  the  program,  means  that  the  value  of  e  is  assigned  to  x.  Such  an 
assignment  also  means  that  no  variable  except  x  is  to  be  altered.  The  concurrent 
assignment,  denoted 

xl,  x2,  ...,  xn  :■  el,  e2,  ...,  en 

means  that  expressions  el,  e2,  ...,  en  are  evaluated  independently,  and  their 


values  assigned  simultaneously  to  xl,  x2 . xn,  respectively.  As  before, 

the  absence  of  a  variable  on  the  left  side  means  that  it  Is  unchanged  by  the 
assignment. 

The  conditional  assignment,  denoted 
Cpl  -*■  A1  |  p2  -*■  A2  |  . . .  |  pn  -*•  An) 

where  pi,  p2,  ...,  pn  are  predicates  and  Al,  A2,  ...,  An  are  assignments 
(simple,  concurrent  or  conditional)  means  that  particular  assignment  Ai 
associated  with  the  first  pi.  If  any,  which  evaluates  true;  otherwise.  If  no 
pi  evaluates  true,  then  the  conditional  assignment  is  undefined. 

An  expression  in  an  assignment  may  contain  a  function  value,  e.g., 
x  :•  max  (x,  abs(y)) 

where  max  and  abs  are  functions.  But  the  function  defined  by  the  assignment 
statement  is  different,  of  course,  from  max  or  abs. 

We  note  that  many  programming  languages  permit  the  possibility  of  so- 
called  side  effects,  which  alter  data  not  mentioned  in  assignment  statements 
or  in  tests.  Side  effects  are  specifically  prohibited  in  our  definition  of 
assignments  and  tests. 

Proper  programs  -  We  define  a  proper  program  to  be  a  program  whose  flow¬ 
chart  has  exactly  one  entry  line,  one  exit  line,  and,  further,  for  every  node 
a  path  from  the  entry  through  that  node  to  the  exit.  For  example. 


HZb-  "HMD* 


are  not  proper  programs 


Program  functions  -  We  define  a  program  function  of  a  proper  program  P, 
denoted  [P],  to  be  the  function  computed  by  all  possible  executions  of  P  which 
start  at  Its  entry  and  terminate  at  Its  exit.  That  is,  a  program  function  [p] 
is  a  set  of  ordered  pairs,  the  first  member  being  a  state  of  the  data  on  entry 
to  P,  the  second  being  the  resulting  state  of  the  data  on  exit.  Note  that 
the  state  of  data  includes  input,  output  files  which  may  be  read  from  or 
written  to  intermittently  during  execution.  Also  note  that  if  a  program  does 
not  terminate  by  reaching  its  exit  line  from  some  initial  data  at  its  entry, 
say  by  looping  indefinitely  or  by  aborting,  no  such  pair  will  be  determined  and 
no  trace  of  this  abnormal  execution  will  be  found  in  its  program  function. 

Proper  programs  are  convenient  units  of  documentation.  Their  program 
functions  abstract  their  entire  effect  on  the  data  known  to  the  program. 

Within  a  program,  any  subprogram  which  is  proper  can  be  also  abstracted  by  its 
program  function,  that  is,  the  effect  of  the  subprogram  can  be  described  by  a 
single  function  node  whose  function  is  the  program  function  of  the  subprogram. 

We  say  two  programs  are  function  equivalent  if  their  program  functions 
are  identical.  For  example,  the  programs 


have  different  flowcharts  but  are  function  equivalent. 

Prime  programs  -  We  define  a  prime  program  to  be  a  proper  program  which 
contains  no  subprogram  which  is  proper  except  for  itself  and  function  nodes. 
For  example. 


8 


*;  s 


if  p  then  f  fi 


while  p  do  f  od 


do  f  until  p  od 


if  p  then  f  else  g  fi 


dol  f  while  p  do 2  y  od 


Larger  primes  will  go  unnamed  here,  although  the  case  statement  of 
Pascal  is  a  sample  of  a  useful  larger  prime.  All  of  the  primes  above  except 
the  last  (dowhiledo)  are  common  to  many  programming  languages.  Prime  programs 
in  text  form  can  be  displayed  with  standard  indentation  to  make  the  subprogram 
structure  and  control  logic  easily  read,  which  we  will  illustrate  for  ZEROIN. 


Y-yr 
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III.  UNDERSTANDING  ZEROIN 

The  prime  program  decomposition  of  ZEROIN  -  Our  first  step  in  under¬ 
standing  ZEROIN  was  to  develop  a  prime  program  decomposition  of  Its  flowchart. 
After  a  little  experimentation,  the  flowchart  for  ZEROIN  was  diagrammed  as 
shown  in  Figure  2.  The  numbers  in  the  nodes  of  the  flowchart  represent 
contiguous  segments  of  the  FORTRAN  program  of  Figure  1,  so  all  lowest  level 
sequence  primes  are  already  identified  and  abstracted. 

The  flowchart  program  of  Figure  2  was  then  reduced,  a  step  at  a  time,  by 
identifying  primes  therein  and  replacing  each  such  prime  by  a  newly  numbered 
function  node,  e.g.,  R. 2.3  names  prime  3  in  reduction  2  of  the  process.  This 
reduction  is  shown  in  Figure  3,  leading  to  a  hierarchy  of  6  levels.  Of  all 
primes  shown  in  Figure  3,  we  note  only  two  which  contain  more  than  one  predi¬ 
cate,  namely,  R.3.1  and  R.S.l,  and  each  of  these  is  easily  modified  into  a 
composite  made  up  of  primes  with  no  more  than  one  predicate.  These  modifica¬ 
tions  are  shown  in  Figure  4.  We  continue  the  reduction  of  these  new  composite 
programs  to  their  prime  decompositions  in  Figure  5.  In  each  of  these  two  cases, 
a  small  segment  of  programs  is  duplicated  to  provide  a  new  composite  which 
clearly  executes  identically  to  the  prime.  Such  a  modification  which  permits 
a  decomposition  into  one  predicate  primes  is  always  possible,  provided  an 
extra  counter  is  used.  In  this  case,  it  was  fortunate  that  no  such  counter 
was  required.  It  was  also  fortunate  that  the  segments  duplicated  were  small; 
otherwise,  a  program  call  in  two  places  to  the  duplicated  segment  might  be  a 
better  strategy. 

A  structured  design  of  ZEROIN  -  Since  a  prime  program  decomposition  of  a 
program  equivalent  to  ZEROIN  has  been  found  with  no  primes  of  more  than  one 
predicate,  we  can  reconstruct  this  program  in  text  form  in  the  following  way: 
The  final  reduced  program  of  ZEROIN  is  given  in  Reduction  6  of  Figure  3,  namely. 


ZEROIN 


Reduction  3 


Figure  3  (3  of  4  pages) 


Reduction  6 


R.6.1  * 


Figure  3  (4  of  4  pages) 


can  be 

modified 

to 


can  be 

modified 

to 


Figure  4 


Reduction  4 
T.4.1 


Figure  5  (2  of  2  pages) 
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Chat  R. 6.1  Is  a  sequence,  repeated  here. 


R. 6.1  • 

I  R. 2.1  | 

T 

rrfn 

1100-102  I 

1  l 

Now  R.2.1  can  be  looked  up,  In  turn,  as: 

K-2-1  ’  , _ 

I  ^  I 


R.1.1 


etc.,  until  all  Intermediate  reductions  have  been  eliminated.  Recall  that 
R.5.1  (and  R. 3.1)  was  further  reduced  in  Figure  5.  When  these  intermediate 
reductions  have  all  been  eliminated,  we  obtain  a  structured  program  in  PDL 
(Process  Design  Language)  for  ZEROIN  shown  in  Figure  6.  Note  there  are  three 
columns  of  statement  numberings.  The  first  column  holds  the  PDL  statement 
number;  the  second  holds  the  FORTRAN  line  numbering  of  Figure  1;  the  third 
holds  the  FORTRAN  statement  numbering  of  Figure  1.  The  FORTRAN  comments  have 
been  kept  intact  in  the  structured  program  and  appear  within  square  brackets 
[,].  From  here  on,  statement  numbers  refer  to  the  PDL  statements  of  Figure  6. 

The  duplication  of  code  introduced  in  Figure  4  can  be  seen  in  PDL  72,  73, 
and  PDL  96-99.  It  should  be  noted,  however,  that  in  PDL  87-91  the  second 
IF  STATEMENT  in  FORTRAN  93  can  be  eliminated  by  use  of  the  if-then-else.  This 
permits  an  execution  time  improvement  to  the  code.  A  second  improvement  can  be 
seen  in  PDL  62-66.  The  use  of  the  absolute  value  function  can  be  eliminated 
and  the  if-then-else  can  be  used  to  transform  the  else  negative  p  into  a 
positive  p  only  in  the  case  where  p  is  negative. 


FORTRAN 

.ine 

Scat 

lefer- 

1 

ace 

! 

ref. 

ZEROIN.  PROGRAM 

i 

1-2 

func  zeroin  (real  ax,  bx,  f.  Col.  integer  Id) 

2 

5 

real  a,  b,  c,  d,  e,  eps,  fa,  fb,  fc. 

3 

Col  1,  xm,  p,  q,  r,  s 

4 

7 

[COMPUTE  EPS,  THE  RELATIVE  MACHINE  PRECISION] 

5 

9 

eps  :•  1.0 

6 

do 

7 

10 

10 

eps  :■  eps/2.0 

8 

11 

col  1  :■  1.0  +  eps 

9 

1 

unCll 

10 

12 

col  1  <  1 

11 

od 

12 

14 

1 

[INITIALIZATION] 

13 

16 

if  ip  -  1  Chen  write  (’THE  INTERVALS  DETERMINED 

14 

18 

a  !■  ax 

15 

19 

b  bx 

16 

20 

fa  :•  f(a) 

17 

21 

fb  f (b) 

18 

23 

[BEGIN  STEP] 

19 

25 

20  1 

c  :■  a 

20 

26 

fc  fa 

21 

27 

d  b-a 

22 

28 

e  :*  d 

23 

dol 

24 

29 

if  ip  *  1  then  write  (b,  c)  fi 

25 

1  30 

if 

26 

31 

abs  (fc)  <  abs  (fb) 

27 

then 

28 

32 

a  :»  b 

29 

33 

b  :■  c 

30 

34 

1 

c  :»  a 

31 

35 

fa  ;«  fb 

32 

36 

fb  fc 

33 

37 

fc  fa 

34 

fi 

35 

39 

[CONVERGENCE  TEST) 

36 

41 

40 

tol  1  :■  2.0  *  eps  *  abs  (b)  +  0.5  *  tol 

37 

42 

xa  :•  .5  *  (c-b) 

38 

to 

while 

39 

abs  (xa)  >  tol  1  and  fb  ^  0 

40 

do2 

41 

{IS  BISECTION  NECESSARY] 

42 

if 

43 

dbi^e)  <  tol  1  or  abs  (fa)  &  abs  (fb) 

44 

83 

70 

then  [BISECTION] 

45 

85 

d  xa 

46 

86 

e  :■  d 

47 

46 

else  (IS  QUADRATIC  INTERPOLATION  POSSIBLE] 

48 

if 

49 

48 

a  i  c 

50 

62 

then  [INVERSE  QUADRATIC  INTERPOLATION] 

fi 


Figure  6.  (1  of  2  pages) 


Lins  Stmt 
Refer-  # 
•nee  ref. 


51 

64 

52 

65 

53 

66 

54 

67 

55 

68 

56 

55 

57 

57 

58 

58 

59 

59 

60 

61 

70 

62 

63 

72 

64 

65 

72 

66 

67 

73 

68 

75 

69 

70 

77 

71 

83 

72 

85 

73 

86 

74 

75 

79 

76 

80 

77 

78 

79 

80 

90 

81 

91 

82 

83 

92 

84 

85 

92 

86 

87 

88 

93 

89 

90 

93 

91 

q  :•  fa/fc 
r  :•  fb/fc 
•  :■  fb/fa 

p  :•  a  *  (2.0  *  jem  *  q  *  (q-r)  -  (b-a)  *  (r-1.0)) 
q  :«  (q-1.0)  *  (r-1.0)  *  (s-1.0) 
alia  [LINEAR  INTERPOLATION] 
a  :•  fb/fa 
p  :•  2.0  *  xm  *  a 
q  :■  1.0  -  a 


fi 

Tajd just  signs] 

if 

p  >  o 
than 

q  -q 
fi 


/*  note  can  be 
/*  .if  p  >  o  then  q  :■ 
/*  else  p  :» 

/*  in  PDL 
/* 


p  abs(p) 

[IS  INTERPOLATION  ACCEPTABLE] 
if 

(2.0  *  p)  2  (3.0  *  xm  *  q  -  abs  (tol  1  *  q)) 
then  [BISECTION] 

d  xm  /*  note  85-86  repeated 

e  :■  d  /*  in  PDL 


else 

a  :»d 


*/ 
■q  */ 

■P  */ 


*/ 

*/ 


d  :•  p/q 
fi 
fi 

[COMPLETE  STEP] 

a  :«  b 


fa  :«  fb 


if 

aba(d)  >  tol  1 
then 

b  :■  b  +  d 
fi 
if 

eba(d)  5  tol  1 
than 

b  :«  b  +  aign  (tol  1,  xm) 
fi 


/* 

note  test  done  twice 

*/ 

/* 

in  FORTRAN 

*/ 

/* 

here  and 

*/ 

f* 

in  line  88 

*/ 

92 

93 

94 

95 

96 

97 

98 

99 
100 


94 


25 

26 

27 

28 


101 

102 

103 

104 

105 


98 

100 

101 

102 


20 


fb  :•  f(b) 
if 

fb  *  (fc/aba  (fc))  >  0.0 
than  [BEGIN  STEP] 
c  :•  a 
fc  :•  fa 
d  :■  b  -  a 
a  :■  d 
fi 
od 

Toone] 
saroin  :•  b 
return 
enuf 


/*  note  25-28  *f 
/*  repeated  */ 
/*  in  PDL  */ 


Figure  6.  (2  of  2  pages) 
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By  construction,  the  PDL  program  of  Figure  6  is  function  equivalent 
to  the  FORTRAN  program  of  Figure  1.  But  the  PDL  program  will  be  simpler 
to  study  and  understand. 

Data  references  in  ZEROIN  -  Our  next  step  in  understanding  ZEROIN  was  to 
develop  a  data  reference  table  for  all  data  Identifiers.  While  straight¬ 
forward  and  mechanical,  there  is  still  much  learning  value  in  carrying  out 
this  step,  in  becoming  familiar  with  the  program  in  the  new  structured  form. 

The  results  are  given  in  Figure  7.  This  familiarization  led  to  the  following 
observations  about  the  data  references  in  ZEROIN  (in  no  particular  order  of 
significance,  but  as  part  of  a  chronological,  intuitive,  discovery  process): 

1.  ax,  bx,  f,  ip,  tol  are  never  set,  as  might  be  expected,  since 
they  are  all  input  parameters  (but  this  check  would  determine 
initialized  data  if  it  existed,  and  also  checks  for  the  presence  of 
side  effects  by  the  program  on  its  parameters  if  passed  by  reference) . 

2.  Zero in  is  never  used,  but  is  returned  as  the  purported  zero  found 

for  f  (since  Zeroln  is  set  to  b  just  before  the  return  of  the  program, 
it  appears  that  b  may  be  a  candidate  for  this  zero  during  execution). 

3.  eps  is  set  by  the  dountil  loop  6-11  at  the  start  of  program  execution, 
then  used  as  a  constant  at  statement  36  from  then  on. 

4.  tol  1  is  used  for  two  different  unrelated  purposes,  namely,  as  a 
temporary  in  the  dountil  look  6-11  which  sets  eps,  then  reset  at 
statement  36  as  part  of  a  convergence  consideration. 

5.  the  function  f  is  called  but  three  times,  $4^16,  17  to  initllize 

fa,  fb,  and  at  92  to  reset  fb  to  f(b)  (m/re  evidence  that  b  is  the 

/ 

candidate  zero  to  be  returned).  / 

/ 

6.  the  Identifiers  a,  c  are  set  to  and  from  b,  and  the  triple  a,  b,  c 
seems  to  be  a  candidate  for  bracketing  the  zero  which  b  (and  zeroln) 
purports  to  approach. 


a 

14,28,80 

16,19,21,30,49,54,96,98 

ax 

14 

b 

15,29,83,90 

17,21,24,28,36,37,54,80,85,90,92,98,103 

bx 

15 

c 

19,30,96 

29,37,4 9 

d 

21,45,72,76,98 

22,46,73,75,83,85,88,99 

e 

22,46,73,75,99 

43 

eps 

5,7 

7,8,36 

f 

16,17,92 

fa 

16,31,81 

20,33,43.51,53,57,97 

fb 

17,32,92 

26,31,39,43,52,53,57,81,94 

fc 

20,33,97 

26,32,51,52,94 

ip 

13,24 

P 

54,58,67 

63,67,70,76 

q 

51,55,59,65 

54,55,65,70,76 

r 

52 

54,55 

3 

53,57 

54,55,58,59 

col 

36 

col  1 

8,36 

10,39,43,70,83,88 

xa 

37 

39,45,54,58,70,72,90 

zero la 

101 

Figure  7. 

7.  Che  identifiers  fa,  fb,  fc  are  evidently  standins  for  f(a),  f(b), 
f(c),  and  serve  to  limit  the  calls  on  function  f  to  a  minimum. 

8.  the  identifiers  p,  q,  r,  s  are  Initialized  and  used  only  in  the 
section  of  the  program  that  the  comments  Indicate  is  concerned 
with  interpolation. 

9.  focusing  on  b,  aside  from  initialization  at  statement  15,  and  as 
part  of  a  general  exchange  among  a,  b,  c  at  statement  28-29,  b  is 
updated  only  in  the  ifthenelse  83-90,  incremented  by  either  d  or 
tol  1. 

10.  d  is  set  to  xm  or  p/q  (as  a  result  of  a  more  complex  bisection  and 
interpolation  process);  xm  is  set  only  at  statement  37  to  the  half 
interval  of  (b,  c)  and  appears  to  give  a  bisection  value  for  b. 

A  function  decomposition  of  2ER0IN  -  The  prime  program  decomposition  and 
the  familiarity  developed  by  the  data  reference  tabulation  and  observations 
suggest  the  identification  of  various  intermediate  prime  or  composite  programs 
in  playing  Important  roles  in  summing  up  a  functional  structure  for  ZEROIN. 

Each  such  intermediate  prime  or  composite  program  computes  values  of  a  function. 
The  inputs  (function  arguments)  of  this  function  are  defined  by  the  initial 
values  of  all  identifiers  which  are  inputs  (function  arguments)  for  statements 
which  make  up  the  intermediate  program.  The  outputs  (function  values)  of 
this  function  are  defined  by  the  final  values  of  all  identifiers  which  are 
outputs  (function  values)  for  statements  which  make  up  the  intermediate  pro¬ 
gram.  Of  course,  further  analysis  may  disclose  that  such  a  function  is 
independent  of  some  inputs,  if,  in  fact,  such  an  identifier  is  always 
initialized  in  the  intermediate  program  before  its  use. 

On  the  basis  of  this  prime  decomposition  and  data  analysis,  we  reformulated 
ZEROIN  of  Figure  6  as  zeroinl,  a  sequence  of  four  intermediate  programs,  as 
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shown  In  Figure  8,  with  function  statements  using  the  form  f.  n-m  where  n,  m 
are  the  boundary  statements  of  the  intermediate  programs  of  ZERO IN  from  Figure  6. 
The  identifier  *outfile  in  the  output  lists  refers  to  the  fact  that  data  is 
being  transferred  to  an  outflle  by  an  intermediate  program.  The  phrase  (x,z,v) 
projection  of  some  function  x,y,z,u,v,w  :■  p,q,r,s,t,u  means  the  new  function 
x,z,v  p,r,t. 

In  the  program  descriptions  which  follow,  all  arithmetic  operations 
are  assumed  to  represent  machine  arithmetic.  However,  we  will  occasionally 
apply  normal  arithmetic  axioms  in  order  to  simplify  expressions.  We  next 
look  at  the  intermediate  programs. 

f . 5-11  -  The  Intermediate  program  which  computes  the  values  of  f.5-11 

t 

is  a  sequence,  namely,  an  initialized  dountil,  l.e. 

5  eps  :■  1.0 

6  do 

7  eps  :*  eps/2.0 

8  tol  1  :■  1.0  +  eps 

9  until 

10  tol  U  1 

11  od 

After  some  chinking,  we  determined  chat  at  PDL  6,  an  invariant  of  the  form 
16  ■  (3  k  >  o  (eps  •  2  S)  A  1  +  eps  >  1 
must  hold,  since  entry  to  PDL  6  must  come  from  PDL  5  or  PDL  10  (and  in  the 
latter  case  tol  1  >  1,  having  just  been  sec  to  1.0  +  eps,  so  1.0  +  eps  >  1). 
Furthermore,  at  PDL  9  the  invariant 

19  -  O  k  >  1  (eps  -  2_k>)  A  col  1  -  1  +  eps 
must  hold,  by  observing  the  effect  of  PDL  7,  8  on  Che  invariant  16  at  PDL  6. 
Therefore,  at  exit  (if  ever)  from  the  segment  PDL  5-11,  we  must  have  the 
condition  19  A  PDL  10,  namely 

(Jk  >  1  (eps  -  2~k) )  A  1  +  2  eps  >  1  A  toll  •  1  +  «ps  d  1 


1  func  zeroin  1  (real  ax,  bx,  f,  tol.  Integer  ip) 

2  real  a,  b,  c,  d,  e,  eps,  fa,  fb,  fc,  p,  q,  r,  s,  tol  1,  xm 

3  Integer  ip 

4  [compute  eps,  the  relative  machine  precision] 

5  eps,  tol  1  f.  5-11 

6  [initialize  data] 

7  a,  b,  c,  d,  e,  fa,  fb,  fc,  *outfile  f.  13-22  (ip,  ax,  bx,  f) 

8  [estimate  b  as  a  zero  of  f] 

9  a,  b,  c,  d,  e,  fa,  fb,  fc,  p,  q,  r,  s,  tol  l,xa,  *outfile  :■ 

f.  23-101  (a,  b,  c,  d,  e,  f,  fa,  fb,  fc,  ip,  p,  q,  r,  s,  tol  1,  xm) 

10  [set  zeroin  for  return,  zeroin  :•  b] 

11  zeroin  5-  f.  103-103(b) 

12  return 

13  cnuf 


Figure  8 
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Thus  we  have 

Lemma  3-11  The  program  function  of  f.5-11  is  the  constaat  function. 

{(0,  (eps,  tol  1))  l  (3  k  *  1  (eps  •  2  ^))  A  1  +  2  eps  >  1  A  tol  1*1+  eps  ^  1} 
Since  tol  1  Is  reassigned  (in  PDL  36)  before  it  is  used  again,  f.5-11  can  be 
thought  of  as  computing  only  eps. 

f .13-22  -  The  intermediate  program  which  computes  the  value  of  f. 13-22 
is  a  sequence  which  can  be  written  directly  as  a  multiple  assignment.  It  Is 
convenient  to  retain  the  single  output  statement  PDL  13,  and  write 
f. 13-22  -  f. 13-13;  f. 14-22 

yielding 

Lemma  13-22  The  (a,b,c,d,e,*outf ile)  projection  of  f. 13-22  is  function 
equivalent  to  the  sequence 

f. 13-13;  f. 14-22 

where  f. 13-13  -  if  ip  -  i  then  write  ('THE  INTERVALS  DETERMINED  BY  ZEROIN  ARE') 
f. 14-22  «  a,b,c,d,e  :•  ax,bx,ax,bx-ax,bx-ax 
f. 23-101  -  The  intermediate  program  which  computes  the  value  of  f. 23-101 
is  a  bit  more  complicated  than  the  previous  program  segments  and  will  be  broken 
down  into  several  subsegments.  We  begin  by  noticing  that  several  of  the  input 
and  output  parameters  may  be  eliminated  from  the  list.  Specifically,  as  noted 
earlier,  p,  q,  r,  and  s  are  local  variables  to  f. 23-101  since  they  are  always 
recalculated  before  they  are  used  in  f. 23-101  and  they  are  not  used  outside  of 
f. 23-101.  The  same  is  true  for  xm  and  tol  1.  fa,  fb,  and  fc  can  be  eliminated 
since  they  are  only  used  to  hold  the  values  of  f(a),  f(b)  and  f(c). 

After  considerable  analysis  and  a  number  of  false  starts  leading  into  a 
great  deal  of  detail,  we  discovered  an  amazing  simplification,  first  as  a  con¬ 
jecture,  then  as  a  more  precise  hypothesis,  and  finally  as  a  verified  result. 

This  simplification  concerned  the  main  body  of  the  iteration  of  zeroin,  namely 
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PDL  41-92,  and  obviated  the  need  to  know  or  check  what  kind  of  interpolation 
strategy  was  used,  step  by  step.  This  discovery  was  that  the  new  estimate 
of  b  always  lay  strictly  within  the  interval  bracketed  by  the  previous  b  and  c. 
That  Is,  PDL  41-92,  among  other  effects,  has  the  (b)  projection 
b  :■  b  +  o(c-b),  for  some  a,  0  <  a  <  1 
so  that  the  new  b  was  a  fraction  a  of  the  distance  from  the  previous  b  to  c. 

With  a  little  more  thought,  it  became  clear  that  the  precise  values  of  d,  e 
could  be  ignored,  their  effects  being  captured  in  the  proper  (but  precisely 
unknown)  value  of  a.  Furthermore,  this  new  indeterminate  (but  bounded)  variable 
a  could  be  used  to  summarize  the  effect  of  d,  e  in  the  larger  program  part 
PDL  23-101,  because  d,  e  are  never  referred  to  subsequently.  Thus,  we  may 
rewrite  f. 23-101  at  this  level  as 

a,  b,  c  *outfile  :■  f. 23-101  (a,  b,  c,  f,  ip) 
and  we  define  it  as  an  Initialized  while  loop. 

Lemma  23-101  The  (a,  b,  c,  *outfile)  projection  of  f. 23-101  is  function 
equivalent  to 


(ip  -  1  •*  write  (b,  c)  |  true  -*•  I); 

[Lemma 

24] 

(  |  f  (c)  |  <  |  f  (b)  |  a,  b,  c  b,  c,  b  [  true  *►  I) 

[Lemma 

25-34] 

while 

f(b)  i  0  A  1  (c-b)/2  |  >  2  eps  |  b  |  +  tol/2 

do 

a,  b,  c  :•  b,  b  +  a(c-b),  c  where  0  <  a  <  1; 

[Lemma 

41-92] 

(f(b)  *  f(3>  0  *►  a,  b,  c  a,  b,  a  |  true  -*■  I); 

[Lemma 

93-100] 

(ip  ■  1  •*  write  (b,  c)  |  true  ♦  1); 

[Leona 

24] 

(  |  f(c)  |  <  |  f(b)  1  a,  b,  c  :*  b,  c,  b  |  true  -►  I) 

[Lemma 

25-34] 

od 

where  1  is  the  identity  mapping. 


I 
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The  structure  of  f. 23-101  corresponds  directly  to  the  structure  of 
POL  23-101  except  for  e  duplication  of  segment  PDL  23-34  in  order  to  convert 
the  dovhiledo  into  a  whiledo.  The  proof  of  the  correctness  of  the  assignments 
of  f. 23-101  is  given  in  separate  lemmas  as  noted  in  the  comments  attached  to 
the  functions  in  Lemma  23-101.  The  while  test  is  obtained  by  direct  substitu¬ 
tion  of  values  for  tol  1  and  xin  defined  in  PDL  36-37  into  the  test  in  PDL  39 
using  eps  as  defined  in  Lemma  5-11. 

Lemma  24  PDL  24  is  equivalent  to  (ip  ■  1  -*■  write  (b,  c)  |  true  -*•  I) 
pf :  By  direct  inspection 

Lemma  25-34  The  (a,  b,  c)  projection  of  the  program  function  of  PDL  25-34 
is  function  equivalent  to 

(  |  f (c)  |  <  |  f (b)  1  -  a,  b,  c  :«  b,  c  b  |  true  -  I) 
pf :  By  direct  inspection  of  PDL  25-34 

Lemma  41-92  The  (a,  b,  c)  projection  of  the  program  function  of  PDL  41-92 
is  function  equivalent  to 

a,  b,  c  :■  b,  b  +  a(c-b),  c  where  0  <  a  <  1 
The  proof  will  be  done  by  examining  the  set  of  relationships  that  must 
hold  among  the  variables  in  PDL  41-92  and  analyzing  the  values  of  p  and  q  only. 
That  is,  it  is  not  necessary  to  have  any  knowledge  of  which  Interpolation  was 
performed  to  be  able  to  show  that  Che  new  b  can  be  defined  by 
b  :•  b  +  a(c-b)  ,  0  <  a  <  1 

We  will  Ignore  the  test  on  PDL  48  since  it  will  be  Immaterial  to  the  lemma 
whether  linear  or  quadratic  interpolation  is  performed.  We  will  examine  only 
the  key  tests  and  assignments  and  do  the  proof  in  two  basic  cases— lncerpo- 
latlon  and  bisection— to  show  that  the  (d)  projection  of  the  program 
function  of  PDL  41-78  is 

d  •  (c-b)  (a)  where  0  <  a  <  1 
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Case  1  Interpolation 

If  Interpolation  is  done,  an  examination  of  Figure  6  shows  that  the 


following  set  of  relations  holds  at  PDL  78: 

•  II.  tol  1  »  2  *  eps  *  abs  (b)  +  .5  *  tol  (PDL  36) 

•  12.  xm  -  (c-b)/2  (PDL  37) 

•  13.  abs  (xm)  >  tol  1  (PDL  39) 

•  14.  p  *  0  (PDL  67) 

*15.  2.  *p<3*xm*q-  abs(tol  1  *  q)  (PDL  70) 

•  16.  d  ■  p  /  q  (PDL  76) 

•  17.  abs(d)  >  tol  1  (PDL  83) 


Now  let ' s  examine  the  set  of  cases  on  p  and  q 
p  >  0  A  q  <  0 

We  have  d  ■  p/q  <  0  (by  hypotheses) , 

p  >  3  xm  +  tol  1  (by  15),  and  col  1  >  0  by(ll) 
q  2  2 

Since  abs(xm)  >  tol  1  (by  13)  and  ^  xm  +  tol  1  <  0  (since  p/q  <  0) 

we  have  xm  <  0  implying  0  >  d  >  p  >  3  xb  >  3  (c-b)  >  (c-b) . 

q  2  T 

Thus  0  >  d  >  (c-b)  yielding  d  ■  a (c-b)  where  0  <  a  <  1 
P  >  0  A  q  >  0 

We  have  d  «  p  >  0  (by  hypotheses), 

q 

£  <  3  xm  -  tol  1  <  3  xm  -  3  (c-b)  <  (c-b)  (by  15,  11,  12) 
q  I  2  4 

implying  0  <  d  <  (c-b).  Thus  d  ■  a(c-b)  where  0  <  a  <  1 
p  >  0  A  q  ■  0 

q  *  0  implies  0  >  2  *  p  (by  15)  and  we  know  p  >  0  (by  hypotheses), 
implying  a  contradiction 
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p  «  0  A  a  *  anything 

abs(p/q)  >  tol  1  (by  16,  17)  and  Col  1  0  (by  II)  implies  p  cannot 

be  0 

p  <  0  A  q  -  anything 

p  >  0  (by  14)  implies  a  contradiction 
Case  2  Bisection 

If  bisection  Is  done,  an  examination  of  Figure  6  shows  that  the  follow¬ 
ing  set  of  relations  holds  at  PDL  78 

Bl.  xm  -  (c-b)/2  (PDL  37) 

B2.  abs(xm)  >  tol  1  (PDL  39) 

B3.  d  -  xm  (PDL  45  or  PDL  72) 

Here  d  ■  xm  (by  B3)  Implies  a  -  4  (by  Bl)  and  thus  d  »  (c-b) (a)  where 
0  <  a  <  1 

PDL  82-91  implies  If  1  d  |  <  tol  1  (l.e.,  if  d  is  too  small)  then 

Increment  b  by  tol  1  with  the  sign  adjusted  appropriately 

l.e.  set  a  »fd  abs(d) 

1  sign  (tol  1,  xm)  otherwise 

But  tol  1  <  abs(xm)  (by  13  and  B2)  -  abs((c-b)/2)  and  the  sign  (tol  1) 
Is  sec  to  the  sign  (xm)  implying 
col  1  ■  ot (c-b)  where  0  <  a  <  1 


>  tol  l'l 
ise  | 


Thus,  In  PDL  82-91  b  is  Incremented  by  d  or  tol  1,  both  of  which 
are  of  the  form  a (c-b)  where  0  <  a  <  1.  Thus  we  have 


b  b  +  a (c-b)  ,  0  <  o  <  1 

and  since  in  PDL  80-81  we  have  a,  fa  :■  b,  fb  we  get  the  statement  of 


Che  Lemma. 


Once  again,  Che  reader  is  reminded  that  the  proof  of  Lemma  41-92  was  done 
by  examining  cases  on  p  and  q  only.  No  knowledge  of  the  actual  interpolations 
was  necessary.  Only  tests  and  key  assignments  were  examined.  Also,  the  pro¬ 
gram  function  was  abstracted  to  only  the  key  variables  a,  b,  c  and  a  represented 
Che  effect  of  all  other  significant  variables. 

Lemma  93-100  The  (a,b,c)  projection  of  PDL  93-100  is  function  equivalent  to 
(f  (b)  *  f  (c)  >  0  ■*  a,  b,  c  :*  a,  b,  a  |  true  -*•  I) 
pf :  By  direct  inspection,  PDL  93-100  is  an  if  then  statement  with  if 
test  equivalent  to  the  condition  shown  above  and  assignments  which  include 
the  assignments  above. 

The  last  function  in  zeroin  1  (from  Figure  8)  is  the  single  statement 
PDL  103  which  can  be  easily  seen  as 
Lemma  103  f.103  is  function  equivalent  to  zeroin  :■  b 

Now  that  each  of  the  pieces  of  zeroin  1  have  been  defined,  the  program 
function  of  zeroin  will  be  given.  First,  let  us  rewrite  zeroinl,  all  in  one 
place,  using  the  appropriate  functions  (Figure  9). 


1  func  zeroinl  (real  ax,  bx,  f,  tol.  Integer  ip) 

2  real  a,  b,  e,  d,  e,  eps,  fa,  fb,  fc,  a 

3  file  *outfile 

4  [compute  eps,  the  relative  machine  precision] 

5  eps  (x  |  (3  k  >/  1  (x  «  2-k))  Al  +  2x>  1  Al  +  eps  i  1} 

6  [initialize  data] 

7  (ip  -  1  -*■  *outf ile  'THE  INTERVALS  DETERMINED  BY  ZEROIN 

ARE'  |  true  *  I)  ; 

8  a,b,c,d,e  ax,bx,ax,bx-ax,bx-ax 

9  [estimate  b  as  a  zero  of  f] 

10  (ip  ■  1  -*■  *outfile  (b,  c)  1  true  I)  ; 

11  (abs(f(c))  <  abs(f(b))  a,  b,  c  :•  b,  c,  b  j  true  -*■  I) 

12  while 

13  f (b)  *  0  A  |  (c-b)/2  j  >  2  eps  j  b  |  +  tol/2 

14  do 

15  a,  b,  c  :■  b,  b  +  a  (c-b) ,  c  where  0  <  a  <  1;* 

16  (f  (b)  *  f(c)  >  0  -*•  a,  b,  c  a,  b,  a  |  true  I)  ; 

17  (ip  *  1  ■*  *outfile(b,  c)  |  true  ■*  I)  ; 

18  (abs(f(c))  <  abs(f  (b))  ♦  a,  b,  c  b,  c,  b  |  true  I) 

19  od 

20  [set  zeroin  for  return,  zeroin  :»  b] 

21  zeroin  :■  b 

22  return 

23  cnuf 

Figure  9  . 


*  a  is  an  indeterminate  based  on  the  current  values  of  a,  b,  c,  d, 
fa,  fb,  fc,  tol  and  eps 
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Theorem  1-105 

func  zeroln  has  program  function  [zeroin]  - 
(ax  •  bx  -*■  root  bx  | 
f(bx)  ■  0  -*■  root  :»  bx  | 
f  (ax)  »  0  -*■  root  :*  ax  | 

f(ax)  *  f(bx)  <0  -*•  root  :■  approx  (f,  ax,  bx,  tol)  | 

true-*-  (V  k  •  1,2, . . .  ,f  (b^)  *  f(ck)  >  0  root  :*  unpredictable! 

3k  >  0  (f(bk)  *  f(ck)  <  0  A  V  j  -  l,2,...k  -1,  f(b  )  *  fCCj)  >  0)- 
root  :■  approx  (f,  bk>  ck»  tol) 

where 

approx  (f,  ax,  bx,  tol)  Is  some  value  In  the  Interval  (ax,  bx)  within 
4  *  eps  *  |  x  |  +  tol  of  some  zero  x  of  the  function  f 

and 

the  sequence  (bl,  cl),  (b2,  c2),  ...  is  defined  so  that  each 
succeeding  interval  is  a  sub-interval  of  the  preceding  interval; 
and  in  the  case  where  abs(d)4tol  1  never  occurs  (bl,  cl}  »  (ax,  bx}, 

{bk+1  c^^}  defines  the  half  interval  of  (bk>  cfc}  including  bfc,  and 

b.  is  chosen  to  minimize  abs(f (b,  . . )) . 

k+1  krrl 

Proof :  The  proof  will  be  carried  out  in  cases,  corresponding  to  Che  conditions 
in  the  rule  given  in  the  Theorem,  The  first  three  cases  follow  directly  by 
inspection  of  zerolnl,  as  special  cases  for  input  values,  which 
bypass  the  while  loop.  X.e.,  if  ax  ■  bx,  then  the  values  of  a,  b,  c  and 


root  can  be  traced 

in 

zeroln 

1  as 

follows: 

a 

b 

c 

root 

zeroln  1.8 

bx 

bx 

bx 

.11 

bx 

bx 

bx 

(condition  13  fails 

since 

c-b  ■ 

0] 

.21 

bx 

bx 

bx 

bx 

l 
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Cases  2  and  3  proceed  in  a  similar  fashion. 

Case  4,  f(ax)  *  f(bx)  <0,  will  be  handled  by  an  analysis  of  Che  vhlledo 
loop  and  les  results  will  apply  co  Che  last  subcase  of  che  last  case  as  well. 

The  flrsc  subcase  of  Che  last  case  arises  when  no  zero  of  f  Is  even  bracketed 
and  zerolnl  runs  a  predictable  course,  as  will  be  shown. 

Case  4;  It  will  be  shown  that  the  entry  condition  f(ax)  *  f(bx)  <0  leads  to 
the  following  condition  at  the  whiletest  of  zerolnl: 

I«(a»cj*bVa<b<cvc<b<a)  A  f(b)  *  f(c)  5  0  A  abs(f(b)  S  abs(f(c)) 
The  proof  is  by  induction.  First  I  holds  on  entry  to  the  whiledo  loop  because 
by  direct  calculation 

after  zerolnl. 8  a  ■  c  A  f(b)  *  f(c)  <0  A  c  i1  b 

after  zeroinl.il  a  ■  c  A  f(b)  *  f(c)  <0  A  abs(f(b>)  abs(f(c))  A  c  i*  b 

Next,  suppose  che  Invariant  I  holds  at  any  iteration  of  the  whiledo  at 
the  whiletest,  and  the  whiletest  evaluates  true,  it  can  be  shown  that  I  is  pre¬ 
served  by  the  three-part  sequence  of  the  do  part.  In  fact,  it  will  appear  that 
the  first  part,  in  seeking  a  better  estimate  of  a  zero  of  f  may  destroy  this 
invariant,  and  the  last  two  parts  do  no  more  Chan  to  restore  the  invariant. 

It  will  be  shown  in  Lemma  15-18  that 

after  zerolnl. 15  (a  <  b  <  c  v  c  <  b  <  a)  a  f(a)  *  f(c)  <0 

after  zerolnl. 16  (a*c^b  va<b<cvc<b<a)A  f(b)  *  f(c)  SO 

after  zerolnl.  18  (a*ci*b  v  a  <  b  <  c  v  c  <  b  <  a)  A  f(b)  *  f(c)  SO  A 

abs(f(b))  S  abs(f (c)) 

which  is  I,  again.  Thus,  I  is  indeed  an  invariant  at  the  whiletest. 

Consider  che  question  of  termination  of  che  whiledo.  In  Lemma  15-18T 

it  will  be  shown  using  c«  and  b»  as  entry  values  to  the  do  part,  chat  for 

some  a,  0<a<l,  after  zerolnl. 18  abs(c-b)  <  abs(c  -  b  )  max  (a,  1-a) . 

oo 


Therefor*,  Che  whiledo  must  finally  terminate  because  the  condition 
f(b)  +  0  A  abs((c-b)/2)  >  2  *  eps  *  abs(b)  +  tol/2 
must  finally  fail,  because  by  the  finiteness  of  machine  precision  abs(c-b) 
will  go  to  zero  if  not  terminated  sooner. 

When  the  whiledo  terminates,  the  invariant  I  must  still  hold.  In  par¬ 
ticular  f(b)  *  f(c)  {  0,  which  combined  with  the  negation  of  the  whiletest  gives 
IT  -  f (b)  *  f(c)  $  0  A(f (b) )  -  0  V  abs((c-b)/2)  *  2  *  eps  *  abs(b)  +  tol/2 
IT  states  that 

1)  a  zero  of  f  is  bracketed  by  the  interval  (o,  c) 

2)  either  the  zero  is  at  b  or  the  zero  is  at  most  |  c-b  |  from  b, 
i.e.,  the  zero  is  within  A  *  eps  *  |  b  |  +  tol  of  b. 

This  is  the  definition  of  approx  (f,  b,  c,  tol). 

Now,  beginning  with  the  interval  (ax,  bx),  every  estimate  of  b  created  at 
zerolnl.15  remains  within  the  interval  (b,c)  current  at  the  time*.  Since  c 
and  b  are  initialized  as  ax  and  bx  at  zerolnl.8,  the  final  estimate  of  b  is 
given  by  approx  (f,  ax,  bx,  tol).  The  assignment  zeroin  :«  b  at  zerolnl.21 
provides  the  value  required  by  case  4. 

Case  5:  part  1.  We  first  show  that  in  this  case  the  condition  a  *  c  will  hold 
at  zerolnl.15  if  f(b)  *  f(c)  >  0.  By  the  hypothesis  of  case  5,  part  1, 
f((b+c)/2)  is  of  the  same  sign  as  f(b)  and  f(c).  Therefore,  the  first  case  of 
zeroinl.16  will  hold  and  the  assignment  c  :■  a  will  be  executed  implying  a  -  c 
when  we  arrive  at  zerolnl.15  from  within  the  loop.  Also,  if  we  reach  zerolnl.15 
from  outside  the  loop  (zeroinl.8-11)  we  also  get  a  »  c. 

We  now  apply  Lemma  15L,  which  states  that  under  the  above  condition  the 
(a,  b,  c)  projection  of  zerolnl.15  is 

this  is  because  f(b)  *  f(c)  $  0  is  part  of  I 
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(f(b) 


((c)  »  0  -  b,  c  t>.  *  *?>«•  “ 


if  abs(c-b)/2»  tol  1 
otherwise 


trn  a,  b,  c  :•  b,  b  +■  a(c-b),  c) 
which  is  a  refinement  of  zeroinl.15. 

Note  that  zeroinl.18  may  exchange  b,c  depending  on  aba(f(b))  and  abs(f(c)). 
Thus,  the  (b,c)  projection  of  the  function  computed  by  zeroinl. 15-18  in  this 
case  is 

k  _  .  J"b  +  (c-b)/2\  .  u  „  .  fb  +  (c-b)/2\ 

i.e. ,  the  new  interval  (b,  c)  is  the  half  interval  of  the  initial  (b,,  c,) 
which  includes  b0  (for  Increments  greater  than  tol  1),  and  the  new  b  is  chosen 
to  minimize  the  value  abs(f(b)).  The  result  of  iterating  this  dopart  is 
unpredictable  unless  more  is  known  about  the  values  of  f.  For  example,  if  the 
values  of  f  in  (ax,  bx)  are  of  one  sign  and  monotone  increasing  or  decreasing, 
then  the  iteration  will  go  to  the  end  point  ax  or  bx  for  which  abs(f)  is 
minimum.  In  general,  the  iteration  will  tend  coward  a  minimum  for  abs(f),  but 
due  to  Che  bisecting  behavior,  no  guarantees  are  possible. 

Case  5:  part  2.  This  covers  the  happy  accident  of  some  intermediate  pair 
b,c  bracketing  an  odd  number  of  zeroes  of  f  by  happening  Into  values  b^,  c^,  such 
that  f(bjj)  *  ffc^)  -  0*  The  tendency  to  move  towards  a  minimum  for  abs(f(b)) 
may  Increase  the  chances  for  such  a  happening,  but  provide  no  guarantee.  Once 
such  a  pair  b^,  c^  is  found,  case  4  applies  and  some  zero  will  be  approximated. 

This  completes  the  proof  of  the  theorem  except  for  the  proofs  of  the 
three  lemmas  used  in  the  proofs  which  follow  directly. 
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Leama  13-18  The  inverlent  I  defined  as 

I  M*  ■  c  i1  b  V  i  <  b  <  e  V  c  <  b  <  «)  A  f (b)  *  f(c)  $  0  A  abs(f(b))  $  aba(f(c)) 
is  preserved  by  the  execution  of  the  loop  body  ZER0IN1. 15-18. 

Proof:  We  use  the  following  abbreviations: 

P  2  abs(f (b))  i*  0  A  abs((c-b)/2)  >  2  *  eps  *  abs(b)  +  tol/2 

1  5  ((c  <  b)  V  (c  >  b))  A  f(b)  *  f(c)  <  0 

i!  2  (a  <  b  <  c  V  c  <  b  <  a)  A  f(a)  *  f(c)  <  0 

I25(a-c»‘bVa<b<cVc<b<a)A  f(b>  *  f(c)  *  0. 

Note  that  P  is  the  loop  predicate.  The  validity  of  the  Lemma  is  an  immediate 
consequence  of  the  following  conditions: 

Cl  :  I  A  P  =3>  I 

o 

C2  :  I  {ZERO INI. 15}  Ii 
0  * 

C3  :  i!  {ZERO INI. 16}  I2 
C4  :  I2  {ZERO INI. 18}  I 

Condition  Cl  is  straightforward.  C2  can  be  seen  by  considering  c  <  b  and 

c  >  b  as  different  input  cases.  Condition  C3  follows  from 

Ij  A  f(b)  *  f(c)  >0  {c  :•  a}  I2  (note  that  setting  c  ■  a  changes  the 

sign  of  f(c>) 

lx  A  f(b)  *  f(c)  «  0  I2 
Similarly,  C4  can  be  inferred  from 


I2  A  abs(f (c))  <  abs(f (b))  {a,  b,  c  :-  b,  c,  b}  I 

I2  A  abs(f  (c)>  >,  abs(f  (b))  X. 
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Lemma  15-181  Given  b,,  c*  on  entry  to  zeroinl. 15-18  then  for  some  o,  0<a<l 
after  zero ini. 15  abs(c-b)  ■  (1-a)  abs(c*-b.) 
after  zeroinl.16  abs(c-b)  <  abs(c.-b.)  max  (a,  1-a) 
after  zeroinl.18  abs(c-b)  <  aba(c.-b.)  max  (a,  1-a) 
proof :  after  zeroinl.15 

abs(c-b)  ■  abs(c,-ba-a(c0-b0)  •  abs(c0-b,) (1-a)  0<a<l 

abs(b-a)  -  abs(b0+a(c.-b.)  -  b.)  -  abs  a(c.-b,)  0<a<l 

after  zeroinl.16 

abs(c-b)  <  max  abs(c.-b,)  (1-a),  abs(c„-b0)a) 

<  abs(c,-b0)  max  (a,  1-a) 
after  zeroinl.18 

abs(c-b)  £  abs(c,-ba)  max  (a, 1-a)  since  b  and  c  are  unchanged  or 

exchanged. 


It  should  be  noted  that  in  the  above  discussion,  zerolnl.17  was  ignored 
because  its  effect  on  the  calculation  of  the  root  and  termination  of  the  loop 
is  irrelevant. 

Ue  have  one  last  lemma  to  prove. 

Lemma  15L  Given  a  ■  c  and  f(a)  *  f(b)  >0  then  zeroinl.15  calculates  the 
new  b  using  the  bisection  method,  i.e., 

b  :■  b  +  C (b-c)/2  if  abs(c-b)  >  tol  l"| 
l^tol  1  otherwise  J 

proof: 

From  POL  43,  either  abs(f(b))  <  abs(f(a))  or  bisection  is 
done  (PDL  45)  with  d  •  xa  ■  (c-b)/2.  then  PDL  82-91  implies 
b  Tb  +  d  -  b  +  (c-b)/2  if  abs(c-b)/2  >  tol  l^J 

^b  +  tol  1  otherwise  J 


Since  by  hypothesis  a  ■  c,  PDL  49  implies  Inverse  quadratic 
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interpolation  is  not  dons  and  linear  interpolation  (PDL  56)  is 
attempted.  Thus 

s  ■  fb/fa  and  0  <  s  <  1  since  £b  *  fa  >  0  and  abs(fb)  <  abs(fa) 

p  ■  (c-b)  *  s,  using  xa  +  (c-b)/2 

q  ■  1-s,  implying  q  >  o  in  PDL  59 

The  proof  will  be  done  by  cases  on  the  relationship  between  b  and  c. 
c  >  b 

c  >  b  Implies  p  >  0  in  PDL  58.  Since  p  >  0  before  PDL  62,  PDL  65 

sets  q  to  -q,  so  q  <  0.  Then  the  test  at  PDL  70  is  true  since 

2  *  p  •  a.  *  s  is  positive, 

3.0  *  xm  *  q  ■  — ■  (c-b)  *  q  is  negative,  and 
abs(tol  1  *  q)  is  positive 
implying  PDL  70  evaluates  to  true 
and  bisection  is  performed  in  PDL  72-73. 

c  <  b 

c  <  b  Implies  p  <  0  in  PDL  58.  Since  p  <  0  before  PDL  62, 

PDL  65  leaves  q  alone  and  PDL  67  sets  p  >  0  Implying  p  ■  (b-c)  *  x. 

Then  the  test  at  PDL  70  is  true  since 
2  *  p  •  2  *  (b-c)  *  s  is  positive, 

3.0  *  xa  *  q  -  4*'  (c-b)  *  q  is  negative,  and 
abs(tol  1  *q)  is  positive 
implying  PDL  70  evaluates  to  true 
end  bisection  is  performed  in  PDL  72-73. 
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IV.  CONCLUSION 

Answering  the  questions  -  We  can  now  answer  Che  questions  originally 
posed  by  Professor  Vandergrafc. 

Question  1: 

If  the  equation  is  linear,  the  program  will  do  a  linear  interpolation 
and  find  the  root  on  one  pass  through  the  loop,  except  in  the  case  where  the 
size  of  the  Interval  (a,  b)  is  smaller  than  tol  1.  Then  it  will  do  a  bisection 
(from  the  test  at  PDL  43).  Note  the  other  potential  condition  where  it  may 
pass  to  PDL  44  for  bisection  is  if  abs(fa)  -  abs(fb)  (from  PDL  19,  26,  and  43). 
However,  in  this  case  bisection  is  an  exact  solution.  The  case  that  the  size 
of  the  interval  is  smaller  than  tol  1  is  unlikely,  but  can  happen. 

Question  2: 

The  theorem  states  that  if  f(a)  and  f(b)  are  both  of  the  same  sign,  we 
will  get  an  answer  that  is  some  point  between  a  and  b  even  though  there  is  no 
root  in  the  interval  (a,  b)  (case  Sa  of  the  Theorem).  If  there  are  an  even 
number  of  roots  in  the  interval  (a,  b)  then  it  is  possible  the  program  will 
happen  upon  one  of  the  roots  and  return  that  root  as  an  answer  (case  5b  of  the 
Theorem).  To  check  for  this  condition,  we  should  put  a  test  right  at  entry 
to  Che  program  between  PDL  3  and  PDL  4  of  the  form. 

if 

f (a)  *  f (b)  >  0 

then 

write  (’F(A)  and  F(B)  ARE  BOTH  OF  THE  SAME  SIGN,  RETURN  B’) 

else 

PDL  4-102 

fi 


Question  3: 

Ic  would  be  easy  to  remove  the  Inverse  quadratic  Interpolation  part 
of  the  code.  We  can  do  this  simply  by  removing  several  PDL  statements,  i.e. , 

POL  47-55.  However,  this  would  not  leave  us  with  the  best  solution  since 
much  of  the  code  surrounding  the  inverse  quadratic  interpolation  could  be 
better  written.  For  example: 

(1)  there  would  be  no  need  to  keep  a,  b,  and  c 

(2)  the  test  in  PDL  70  could  be  removed  if  we  checked  in  the  loop 
that  f(a)  *  f(b)  was  always  greater  than  zero,  since  bisection 
and  linear  interpolation  would  never  take  us  out ‘of  the  interval. 

Cleaning  up  the  algorithm  would  probably  require  a  substantial  transformation. 
Question  4: 

Zero in  will  find  a  triple  root.  It  will  not  inform  the  user  that  it  is 
a  triple  root,  but  will  return  it  as  a  root  because  once  it  has  a  root 
surrounded  by  two  points  such  that  f(a)  and  f(b)  are  of  opposite  signs,  it 
will  find  that  root  (case  4  of  the  Theorem). 


Program  history  -  Since  most  programs  seen  by  practicing  programmers 
do  not  have  a  history  in  the  literature,  we  did  not  research  the  history  of 
ZEROIN  until  we  had  completed  our  experiment.  The  complexity  of  the  program  is 
partially  due  to  the  fact  that  it  was  modified  over  a  period  of  time  by  differ¬ 
ent  authors,  each  modification  making  ie  more  efficient,  effective  or  robust. 
The  code  is  based  on  the  secant  method  (Ortega  and  Reinboldt).  The  idea  of 
combining  it  with  bisection  had  been  suggested  by  several  people.  The  first 
careful  analysis  seems  to  have  been  by  T.  J.  Dekker  (Dekker). 

R.  F.  Brent  (Brent)  added  to  Dekker's  algorithm  the  Inverse  quadratic  inter¬ 
polation  option,  and  changed  some  of  the  convergence  tests.  The  Brent  book 
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contains  an  ALGOL  60  program.  The  FORTRAN  program  of  Figure  1  Is  found  In 
(Forsythe,  Malcolm  &  Moler)  and  Is  a  direct  translation  of  Brent's  algorithm, 
with  the  addition  of  a  few  lines  that  compute  the  machine-rounding  error. 

He  understand  that  ZEROIN  is  a  significant  and  actively  used  program  for  cal¬ 
culating  the  roots  of  a  function  in  a  specific  interval  to  a  given  tolerance. 

Understanding  and  documenting  -  As  it  turns  out,  we  were  able  to  answer 
the  questions  posed  and  discover  the  program  function  of  ZEROIN.  The  techniques 
used  included  function  specification,  the  discovery  of  loop  invariants,  case 
analysis,  and  the  use  of  a  bounded  indeterminate  auxiliary  variable.  The 
discovery  process  used  by  the  authors  was  not  as  direct  as  it  appears  in  the 
paper.  There  were  several  side  trips  which  included  proving  the  correctness 
of  the  inverse  quadratic  interpolation  (an  interesting  result  but  not  relevant 
to  the  final  abstraction  or  the  questions  posed) . 

There  are  some  implications  that  the  algorithm  of  the  program  was  over- 
designed  to  be  correct  and  that  the  tests  may  be  more  limiting  than  necessary. 
This  made  the  program  easier  to  prove  correct,  however. 

Ue  believe  this  experience  shows  that  Che  areas  of  program  specification 
and  program  correctness  have  advanced  enough  to  make  them  useful  in  understanding 
and  documenting  existing  programs,  an  extremely  important  application  today. 

In  our  case,  we  are  convinced  that  without  the  focus  of  searching  for  a  correct¬ 
ness  proof  relating  the  specification  to  the  program,  we  would  have  learned  a 
great  deal,  but  would  have  been  unable  to  record  very  much  of  what  we  learned 
for  others. 

Hamming  pointed  out  chat  mathematicians  and  scientists  stand  on  each  other's 
shoulders,  but  programmers  stand  on  each  other's  toes.  We  believe  chat  will 
continue  to  be  true  until  programmers  deal  with  programs  as  mathematical  objects, 
as  unlikely  as  they  may  seem  to  be  in  real  life,  as  we  have  cried  to  do  here. 
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